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Abstract 

The GHZ distribution, named after Greenberger, Horne and Zeilinger 
(1989) plays an important role in quantum information theory. We dis- 
cuss random variate generation for random vectors with the distribution 
under the RAM and bit models, and in the context of communication, and 
show that all can be achieved in expected bit complexity and expected 
communication complexity 0(n log n), where n is the dimension of the 
hypercube. 
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The GHZ distribution 



The GHZ distribution (Greenberger, Horne and Zeilinger, 1989) is defined on the 
hypercube {0, 1}". It has n parameters 9i, . . . 6'„, which take values in [0, 2tt), 
and represent angles. The probability mass oi x — (xi,...,x„) G {0,1}" is 
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given by 
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Properties of the GHZ distribution are desribed in Gravel's thesis (2011). Sim- 
ulation attempts are reported in Bancal, Branciard and Gisin (2010), Brassard, 
Cleve and Tapp (1999), Regev and Toner (2007), Steiner (2000) and Toner and 
Bacon (2003). 

The purpose of this note is to describe an algorithm for generating a random 
vector with this distribution in expected bit complexity and bit communication 
complexity 0(n log n), where bit complexity charges one unit of cost for each 
bit manipulated, even in the input parameters 6j,l < j < n. That is, the 
simulation proceeds using only a random but finite number of bits of the input. 
More specifically, we discuss random variate generation for p, both for the RAM 
model (in which real numbers can be stored and manipulated at unit cost) and in 
the bit model (in which every bit that is touched, moved, generated, or otherwise 
manipulated is charged). We also address the situation in which one individual 
knows the parameters 6j (and thus aj and /3j), and another individual needs 
to generate the random vector. That second individual can request bits in the 
binary decompositions of the parameters as needed, and every bit requested 
costs one unit. In all these models, a random vector with the GHZ distribution 
can be generated in O(nlogn) expected time. 

Properties of the GHZ distribution. 

For later perusal, it is easy to establish by induction on n that maxa,p(a;) < 1/2. 
It is also useful to rewrite p{x). To this end, define 



and recall that aj + 0j = 1. Let A C {1, . . . , n} be the collection of indices j for 
which Xj = 1, and let B be the remaining indices. Then p{x) can conveniently 
be rewritten as follows: 



aj = cos 9j , Pj = sin 9j , 




This is a decomposition into three terms. Note that 
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defines a probability distribution over {0, 1}". Indeed, it is easily decomposed 
as 

n 

=n(^^(i -^^H)- 

That is, if £,j is a Bernoulli (1 — a^) random variable, then the random vector 
with independent components given by 

has distribution g. Similarly, if we set 

n 

q'{x) = n n = n (^■^■"? + - - ' 

then the random vector with independent components given by 
has distribution q' . 

Random variate generation in the RAM model 

Note the inequality 

p{x) < q{x) + q'{x). 

Since + = 2, we can apply von Neumann's rejection method (see, 

e.g., Devroye, 1986): 

Algorithm 1 Generator for the GHZ distribution 
1: repeat 

2: Generate U uniformly on [0, 1] 

3; Generate independent Bernoulli random variables , . . . , with param- 
eters a\, . . . , a^. 
4: Generate a random bit S. 
5; if 5 = 1 then 
6: set X (Ci,...,^„) 

7: else 

8: SetX^ (l-^i,...,l-^„) 

9: end if 

10: until U{q{X) + q'{X))<p{X). 



The expected number of iterations before exit is precisely two. Within an 
iteration, under the RAM model, all work is 0{n). We assume, of course, 
the availability of a uniform [0, 1] random variate generator, and a random bit 
source. Therefore, the expected time is 0{n). Furthermore, the time taken by 
the algorithm, divided by n is stochastically smaller than a geometric random 
variable with constant mean, so its tails are exponentially decreasing. 
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Random variate generation in the bit model 



Assume next that only a random bit generator is available and that the param- 
eters are only accessible as follows: we can ask for a number of bits of each 
parameter, but will be charged one unit of cost per bit that is revealed. We can 
modify the algorithm of the previous section so that if N is the total number of 
bits (of the parameters) revealed before halting, E{A^} = 0(nlogr?,). 

To warm up for this, consider the problem of generating a Bernoulli ran- 
dom variable with parameter p. If U — O.U1U2 ... is the binary expansion of a 
uniform [0. 1] random variable U, i.e., Ui, U2, ... is our source of unbiased inde- 
pendent random bits, and if p — 0.piP2 ... is the binary expansion of p, then we 
keep comparing bits Ui and pi for i = 1,2,... until for the first time Ui ^ Pi. 
Then, if C/^ = < pi = 1, we return X = 1, and if Ui = 1 > pi ~ 0, we return 
X = 0. It is clear that X — 1 ii and only if C/ < p. Therefore, X is Bernoulli 
(p). The expected number of bits revealed in p is precisely two. The expected 
number of bits needed from our random bit source also is two. Therefore, one 
would assume that the expected number of bits needed from af , . . . , in one 
given iteration for generating the Bernoulli random variables is 2n. However, 
the knowledge of k bits of Oi does not translate into a knowledge of k bits for 
af , so we need a minor modification. To be specific, let p = for a fixed i. 
Let us say that we know the integer part of 9i (which consists of at most three 
bits, since < < 2tt) and k bits of the binary expansion of the factional 
part of 9i, then we know the value of p up to an error not exceeding 1/2'^ by 
virtue of the fact that the derivative of af with respect to 6i is bounded by 
one. Let p{k) be that estimate of p. Then, setting U{k) = O.Ui . . . Uk, so that 
U{k) <U < U{k) -fl/2^ we know that U <p\i U{k) < p{k) - 2/2'= and that 
U > p ii U{k) > p{k) + 2/2*^. Thus the following algorithm can be employed to 
check if [/ < p, and thus, to generate a Bernoulli (p) random variable: 



Algorithm 2 Generator for the GHZ distribution 

1: Set fc ^ 1 

2; repeat forever 

3: Compute U (k) = O.Ui ...Uk and p{k) 

4: if C/(fc) -p(fc) < then 

5; return Y = 1 

6: else if U{k) - p{k) > ^ then 

7: return Y = 

8; else 

9: fc fc + 1 

10: end if 

11: end repeat 



The Y generated by this algorithm is Bernoulli (p). The number of iterations, 
and thus the number of bits requested from 9j (minus those for the integer part), 
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is a random variable, K. We have 



P{i4: > A;} = P <^ \U{k) 



< P 



\u-p\<^ 



< -r- 

- 2k 



Therefore, 



E{i^} = ^ P{i^ > fc} < 5^ min 1, ^ = 4. 

fe=0 fe=o ^ ^ 



Thus, the expected number of iterations is bounded by 4, and the expected 
number of transmitted bits is bunded by three mkore, or 7, for each of the 
required BernoulU random variables. The random vector of Bernoulli random 
variables requires less than 7n expected number of bits from the parameters. 

In the algorithm of the previous section, in each iteration we generate a pair 
{U, X). We explained above how X can be generated using at most An random 
bits in expected value and 7n expected number of bits from the parameters. 
However, U cannot be generated since it is a real number. Note however that 
we do not really need U, but we merely need to generate a Bernoulli random 
variable Y where 

P{Y = 1} = P {U{q{X) + q'{X)) < p{X)} . 

For this, we adapt the method of the previous paragraph. Wc call the fc-th 
approximant of U, q{X) + q'{X) and p{X) the approximant obtained by trun- 
cating U to the first k bits in its binary expansion, and for the functions of 
X, the approximant obtained by calculating the function exactly given that the 
input parameters of the function are truncated to a + A; fractional bits, where 
a > is a given constant. We always send the integer parts of all parameters 9i 
for a fixed total of not more than 3n bits. Wc require that the error in the fc-th 
approximant is at most 1/2^ . We will see that for our functions, we can take 

a = 2 + [log2 n] . 

Call the approximant of U{q{X) + q'{X)) Lk [L for left) and that of p{X) Rk 
{R for right). Then using 6 for possibly different real numbers at most one in 
absolute value, 



\Lk-U{q{X) + q'{X))\ 



< 



< 



q{X) + q'{X) 



e 



2k+c 



U{q{X) + q'{X)) 



U 

2k 

3 

2k- 



2k 



1 



Similarly, 



\Rk-p{X)\<^. 
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Thus, we know that y = 1 if 

3 1 

Also, y = if 

3 1 

^k - ^ > Rk + ^■ 

Thus, the following algorithm can be used to compute the value of Y: 
Algorithm 3 Generator for the GHZ distribution 



Given is X e {0, 1}" 
Set fc ^ 1 
repeat forever 

Compute Lk and Rk 
if Lk-Rk < then 

return Y = 1 
else if Lk — Rk > then 

return Y = 
else 

k + 1 

end if 
end repeat 



Since Lk U {q{X) + q'{X)) and Rk — >■ p{X) as fc t oo, this algorithm halts 
with probability one. If it has not halted after having processed Lk and Rk, 
then we know that 

\Uiq{X) + q'{X)) - p{X)\ < \Lk -Rk\ + ^<^. 

But since U is uniformly distributed on [0, 1], the probability of this is thus not 
more than 



2^^iqiX)+q'iX))- 

Let K be the value of k upon exit from the loop in the algorithm. We have 



I>{K > k\X} < -j: 



2'=(g(X)+g'(X))- 
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Thus, 

OO 

E{K\X} = J2P{K>k\X} 



k=0 



OO / 

< ^minfl, 

k=0 ^ 



< 



2>'{q{X) + q'{X)) 

E 1+ E 



r , . r r .2^{q{X)+q'{X)) 

fe<log2( ,(x)W(X) J fc>'°g2( ,(x)wm j 



18 + log2 



1 



MX) + q'iX) 
Now, we uncondition, and obtain 



E{K} < 18+ y ^M±s:MiogJ--±--) 

< 18 + l0g2 9 X r N 

(by Jensen's inequality) 
= 18 + (n - 1) 
= n+17. 



It is unimportant for what follows, but worth noting nevertheless, that the 
intermediate bound on E{ii'} involves the binary entropy of the distribution 

(9 + g')/2. 

By Wald's identity, the expected number of iterations of the algorithm for 
computing Y is 2E{/\ } and thus 0(n). The expected number of bits transmitted 
or bits manipulated is not more than 

3n + (a + 1) X 0{n) = 0{n logj n). 

We finish the analysis by computing a suitable value for a. First note that the 

partial derivative of p with respect to any 6i is not more than one. Similarly, 
the partial derivative of q + g' with respect to any 6i is not more than four. By 
Taylor's multivariate series expansion, this implies that the change in q+q' value 
due to the truncation to t bits on the part of all parameters (which changes all 
input parameters by at most 1/2^) is at most This is less than 1/2^^ if 

we set 

i = 2+\k + log2 n] . 
This establishes our choice for a. 
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The general GHZ distribution 



The general GHZ distribution on {0, 1}" is of the form 

p{x) = cos^(7) • p'{x) + sin^(7) • p"{x), x e {0, 1}", 

where 

, / ri n \ 

P'ix) - ( n - + (-1)" n sin(e,- - xjn/2)\ 



and 



p"i^) = J I n - ^^•'^z^) - (-1)" n - 



kj=i i=i 



are both like the distribution we defined in the introduction, and 7 is a sum of 3n 
input parameters that are all angles taking values in [0, 2n), and computation is 
mod 27r. For details, see Gravel (2011). It is easy to see that a random variable 
X with distribution p can be obtained as follows: first generate a Bernoulli 
(cos^ 7) random variable Y. If y = 1, let X have distribution p', and if F = 0, 
let X have distribution p" . The only novelty now is the generation of Y. We 
can argue as in the previous sections, and leave it to the reader to check that 
Y can be generated using an expected number of bits (from the 3n parameters) 
that is 0(T; log n). As argued above, a random vector from p' or p" requires 
only O(nlogn) expected bits from the parameters. 



Discussion 

Knuth and Yao (1976) were the first to study the bit complexity for generating 
random integers with given distribution p{x) assuming a perfect Bernoulli source 
of bits. They showed that any algorithm must use an expected number of bits 
at least equal to the entropy ^^p{x) log2(l/p(a;)), and that the best algorithm 
does not use more than three extra bits above the lower bound. The GHZ 
distribution has binary entropy between and n. This provides the theoretical 
backdrop. Our problem, by virtue of the set-up in terms of parameters, whose 
bits must also sparingly be revealed, is slightly different, but related. 

For further results on the bit model in random variate generation, see chapter 
XIV of Devroye (1986) and the recent thesis of Lumbroso (2012). 
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K*^ '. Introduction 



Classical simulation of an entangled quantum state can be thought as the problem of sampling 
the joint discrete probability distribution of the outcomes obtained by people sharing an entangled 
quantum state on which everybody applies locally a von Neumann measurement on his shared part. 
Each person does not know the description of any other measurements belonging to the other parties. 
Each person after having applied his measurement on the subsystem of the state that he shares 
with the others gets a classical outcome. The joint distribution of the outcomes of every parties 
follows a discrete distribution. Simulating entanglement means simulating the discrete probability 
distribution of the quantum state, and, more specifically, simulating with shared or private random 
resources as well as classical communication resources (a classical channel). What does it mean? 
Intuitively, suppose that an entangled quantum state is shared among n persons. A state is entangled 
if it cannot be factored partially or completely into tensor products. In this paper, we are interested 
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in the state l^*) = :^|0") + called the GHZ state in honor of Greenberger, Horne et Zeilinger 

(1989). Each person owns the description of a von Neumann measurement that he can apply to 
his subsystem. Call these n measurements Mj for j € {1, . . . , n}. It is important in what follows 
that the i^^ person does not know Mj for i ^ j. n persons applying jointly and independently 
their measurements mean they apply &j^iMj to l^*). By applying these measurements, they get 
random outputs \hj). The I''* gets |6i), the 2"*^ gets etc. Quantum mechanics asserts these 
random results are eigenvalues of the Mj. The joint distribution of the these results is what we are 
trying to simulate classically. Each local von Neumann measurement is represented by a point on 
the 3-dimensional sphere (called sometimes the Bloch sphere) . A point on the sphere is represented 
by (0, </>) G [0, 27r) x [0, vr). The joint probability distribution is a convex combination of two 
distributions. The coefficients of the convex combination are related to the 9j rcprcnsenting the 
equatorial parts of the Mj and the two distributions are related to the (j)j representing the real parts 
of M,-. 

Equivalcntly, one can think of the n parties transforming |^') by applying locally a unitary 
transormation independently from each other and then measuring in the so-called computational 
basis in quantum computing. By doing this, the basis reprensenting the eigenvectors is then trans- 
formed into the computational basis used in quantum computing. The eigenvectors are then \xj) 
with Xj G {0, 1}. Hence, the n persons measure in the standard basis of computation once they 
have applied U = to the state |^). The axioms of quantum mechanics assert that the joint 

discrete distribution obtained is given by p{x) = |(x|J7|^)p with x = xi ■ ■ ■ Xn € {0,1}"; p is a 
distribution on the hypercube. Properties and the general derivation of the structure of the discrete 
probability distribution p was done in Gravel's thesis (2011). 

Being able of simulating classically p was the subject of many articles. In fact, when saying 
"simulating GHZ" or "simulating entanglement", it means simulating classically p. How ? With 
shared random variables (source of unbiased independently identically distributed bits) and clas- 
sical communication. For the concepts related to communication complexity like the definitions 
of worst-case or average case complexity and shared random variables, sec Kushilevitz and Nisan 
(1997). Originally, there were Einstein, Podolsky and Rosen (1935) who believed they showed that 
quantum mechanics (of Copenhagen) is not complete because it is impossible to reproduce quantum 
correlations with hidden variables (shared random variables) obtained from |^) when n = 2. These 
correlations are the moments of the distribution of p. In Bell (1964), it was shown that it is impossi- 
ble to reproduce these correlations with only shared random variables. Maudlin (1992) proposed to 
add classical communication resource (classical channel) to reproduce those correlations. Maudlin 
made that suggestion without having the knowledge of the theory of complexity used by computer 
scientists. From that moment, the theory of communication complexity became important in order 
to try reproducing p. A sequence of papers were then published on the simulation of p with n = 2 
and n = 3 with average and worst case complexity with different types of measurements, see Bras- 
sard, Cleve and Tapp (1999), Steiner (2000), Toner and Bacon (2003), Regev and Toner (2007), and 
Bancal, Branciard and Gisin (2010). 

Simulating exactly a discrete probability distribution is always possible even if only a source 
unbiased identically independently distributed bits is available. It is important to note that any 
probability distribution can be exactly simulated if one had access to a source of identically inde- 
pendently continuous uniform random variables on the interval (0, 1). Knuth and Yao (1976) were 
the first to study the exact simulation of probability distributions assuming only a source unbiased 
identically independently distributed bits, and their model was termed the random bit model. In 
the random bit model, one can stored and carry algebraic operations on real numbers and, the only 
random source of information he has access, are unbiased identically independently distributed bits. 
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Knuth and Yao also proved that one can simulate exactly any discrete probability distribution using 
a number of random bits that is between the entropy and the entropy plus three bits. We study 
in the random bit model how sampling the GHZ distribution when the parameters are spread out 
among different persons. Each person knows its own parameter completely since he can store real 
numbers, does not know parameters of the other persons and can only learn the knowledge of other 
people's parameters partially by asking a finite number of bits from their binary representations. 



The GHZ distribution 

The GHZ distribution is defined on the hypercube {0, 1}". It has 4n parameters that are n quadru- 
ples {{iOi, (pi, ipi, Ti)}^^i with (cjj, (pi, ipi,T,i) € [0, 27r)^ x {0, 1} for each i G {1, . . . , n}. The probability 
mass of X = (xi, . . . , Xn) G {0, 1}" is given by 



where 



p{x) = cos'^{^)pi{x)+sm'^{^)p2{x), 



1 2 1 2 

-{ai{x) + a2{x)) , p2{x) = -{&i{x) - a2{x)) , 



n 



ai i^) = n - f ^j) ' a2 (x) = JJ - sin {uj - f Xj) . 

Briefly, the products ai and a2 are related to the amplitudes of the entries of the unitary 
transformations. Squares of ai and a2 are themselve discrete probability distributions, p is a convex 

combination pi and p2- The coefficients cos^(^^^) and sin^(^^^) are related to the phases of the 
entries of the unitary transformations. Properties of the GHZ distribution are desribed in Gravel's 
thesis (2011). 

In this paper, we shall concentrate on pi. Why ? Simply because for simulating p, we need to 
simulate a Bernoulli distribution with parameter cos^(2^) in order to branch either on pi or p2 
and it will be clear at the end how to readapt the argument given below in order to simulate such a 
Bernoulli distribution in expected time 0(ri log(ra)). Also, simulating p2 is like pi with one person 
adding tt to its w-parameter. 

From now on, p shall denote pi just before. The purpose of this paper is to describe an algorithm, 
that can be easily turned into a communication protocol, for generating a random vector with 
distribution p in expected bit complexity and bit communication complexity 0(n log n), where 
bit complexity charges one unit of cost for each bit manipulated, even in the input parameters 
•^jj 1 ^ j ^ n. That is, the simulation proceeds using only a random but finite number of bits of the 
input. The expected bit complexity refers to the number of random bits needed for sampling. The 
expected number of communicated bits refer to the number of bits from the fixed parametrs that 
must be communicated to the sampler. More specifically, we discuss random variate generation for 
p, both for the RAM model (in which real numbers can be stored and manipulated at unit cost) and 
in the bit model (in which every bit that is touched, moved, generated, or otherwise manipulated is 
charged) . We also address the situation in which one individual knows the parameters ujj (and thus 
aj and Pj), and another individual needs to generate the random vector. That second individual 
can request bits in the binary decompositions of the parameters as needed, and every bit requested 
costs one unit. In all these models, a random vector with the GHZ distribution can be generated 
in 0(n log n) expected time. 
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Properties of the GHZ distribution. 

It is useful to rewrite p{x). To this end, define 

aj = cos LOj , Pj = sin Uj , 

and recall that a'j + /3j = 1. Let A C {1, . . . , n} be the collection of indices j for which xj = 1, and 
let B be the remaining indices. Then p{x) can conveniently be rewritten as follows: 

\jeA jeB jeA jeB J 



\jeA j€B j&A jeB 3=1 

This is a decomposition into three terms. Note that 

jeA jeB 

defines a probability distribution over {0, 1}". Indeed, it is easily decomposed as 

n 

= n ~ + (1 ~ ■ 

That is, if Xj is a Bernoulli (1 — a|) random variable, then the random vector with independent 
components given by 

(Xi, . . . ,Xn) 

has distribution qi. Similarly, if we set 

n 

Q2{x) = n «j n ^1 = n (^^^^^ + - ^j)^^ - ^i)) = ^^^^^^ 

jeA jeB j=i 
then the random vector with independent components given by 

has distribution q2. 

Random variate generation in the RAM model 

Note the inequality 

Pix) < qi{x) + qiix). 

Since + ^^{x)) = 2, we can apply von Neumann's rejection method (see, e.g., Devroye, 

1986): [algorithm next page] 



4 



Algorithm 1 Generator for the GHZ distribution 



1 

1 


1 CLICdL 
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Generate a random bit S. 




5 


if S" = 1 then 




6 


set X ^ {Xi,...,Xn) 




7 


else 




8 


set X ^{1-Xi,...,l-Xn) 




9 


end if 




10 


until U{qi{X)+q2{X)) <p{X). 





The expected number of iterations before exit is precisely two. Within an iteration, under the 
RAM model, all work is 0{n). We assume, of course, the availability of a uniform [0,1] random 
variate generator, and a random bit source. Therefore, the expected time is 0{n). Fiuthcrmorc, 
the time taken by the algorithm, divided by n is stochastically smaller than a geometric random 
variable with constant mean, so its tails are exponentially decreasing. 

We shall modify and adapt algorithm to get rid of the use of the continuous uniform U (and 
hence its generation) that is not allowed in the random bit model. Once modified, it will be clear 
how to transform the modified algorithm into a communication protocol. 

Random variate generation in the bit model 

Assume next that only a random bit generator is available and that the parameters are only acces- 
sible as follows: we can ask for a number of bits of each parameter (more precisely a number of bits 
of trigonometric functions of each parameter not the parameter themselve as we will see), but will 
be charged one unit of cost per bit that is revealed. 

To warm up for this, consider the problem of generating a Bernoulli random variable with 
parameter p {not the p use for the discrete distribution). If U = O.U1U2 ... is the binary expansion 
of a uniform (0, 1) random variable U, i.e., Ui,U2, ■ ■ ■ is our source of unbiased independent random 
bits, and if p = 0.piP2 ... is the binary expansion of p, then we keep comparing bits Ui and pi for 
i = 1,2,... until for the first time Ui 7^ pi. Then, if Ui = < pi = 1, wc return X = 1, and if 
Ui = 1 > Pi = 0, we return X = 0. It is clear that X = 1 if and only if U < p. Therefore, X is 
Bernoulli (p) . The expected number of bits revealed in p is precisely two. The expected number of 
bits needed from our random bit source also is two if one had to communicate p to someone else 
having access to Ui,U2 ■ ■ ■ 

In algorithm 1 of the previous section, in each iteration we generate a pair {U,X). However, U 
cannot be generated since it is a real number. Note however that we do not really need U, but we 
merely need to generate a Bernoulli random variable Y where 

P{y = 1} = P {U{qi{X) + q2{X)) < p{X)} . 

The event y = 1 is equivalent to accepting in the rejection algorithm 1. For this, we adapt the 
method of the previous paragraph. We call the A;-th approximant of U, qi{X) + 92 (A) and p{X) 
the approximant obtained by truncating U to the first k bits in its binary expansion, and for 

the functions of A, the approximant obtained by calculating the function exactly given that the 
products ai(A) and a2(A) involved in the functions qi, q2 and p are truncated to /c + 2 bits. It 
is clear that the knowledge of the k + 2 bits of the binary expansion of ai(A) and a2(A) implies 
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the knowledge of k bits for qi{X) + q2iX) and p{X). The parameters are technically the angles ojj. 
What really matters to be sent during the communication are the bits from the binary representions 
o/cos(wj) and sm{uij) and not the bits from the uoj 's themselve. For all j = 1, . . . , n, the j^^ person 
can compute locally cos{u}j) and sm{ujj). We require that the error in the A;-th approximant is at 
most 1/2^. 

Call the approximant of U{qi{X) + q2{X)) Lk [L for left) and that of p{X) {R for right). 
Then using 9 for possibly different real numbers at most one in absolute value, 



\Lk-U{qr{X) + q2{X))\ = 



U + 



2k 



qi{X) + q2{X) + ^]- U{qi{X) + q2{X)) 



U q^{X) + q2{X) 1 
3 

- 2^ 



Similarly, 



Thus, we know that y = 1 if 



Also, y = if 



\Rk-p{X)\ < 



2k- 



Lk + ^<Rk-^- 



- ^ > Rk + ^■ 



Thus, the following algorithm can be used to compute the value of Y that depends on the vector 
X and the a^-'s and /3j's. 



Algorithm 2 Generator for the stopping condition (alleviating the need of a continuous uniform 

random variable) 



1 


Given is X G {0, 1}" 


2 


Set A; 1 


3 


repeat forever 


4 


Compute Lk and R^ 


5 


if Lk- Rk < then 


6 


return Y = 1 


7 


else if Lk — Rk> then 


8 


return y = 


9 


else 


10 




11 


end if 


12 


end repeat 



Since L^ U{qi{X) + q2{X)) and Rk p{X) as k'\ oo, this algorithm halts with probability 
one. If it has not halted after having processed L^ and Rk, then we know that 

\U{qi{X) + q2{X)) -p{X)\ < \Lk -Rk\ + ^<^. 

But since U is uniformly distributed on [0, 1], the probability of this is thus not more than 

8 

2>^{qi{X)+q2{X))- 
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Let K be the value of k upon exit from the loop in the algorithm which is the number of bits needed 
from the binary representations U{qi{X) + q2{X) and p{X) in order to reject or accept. We have 

V{K > k\X} < 



2''{qi{X)+q2{X))' 

Thus by denoting /cq = 3 + logg 



1 

qi(X)+q2(X) 



E{K\X} = Y.P{K>k\X} 

< (A;o + l) + l 
= 5 + log2 



qiiX)+q2iX)J ■ 



Now, we uncondition, and obtain 



xe{o,i}" 



qi{x) + q2{x) 



UwD» ^ «l(x) + 92(x) 

(by Jensen's inequality) 
= 5 + (n - 1) 
= n + 4. 



It is unimportant for what follows, but worth noting nevertheless, that the intermediate bound on 
E{i^} involves the binary entropy of the distribution (gi + g2)/2. Also, it might be interesting 
to give some cases for which the knowledge of A; bits is not sufficient (overlapping cases) are for 
examples : 

p{X) <Lk<Rk< U{qi{X) + q2{X)) (1) 
p{X) < Lfe < U{qi{X) + q2{X)) < Ru (2) 
Lk<V<U{qi{X) + q2{X))<Rk (3) 

We therefore need to know n + 4 bits on average from U{qi{X) + q2{X) and p{X) to reject or 
accept X. Computing the left and right approximants and i?^ is therefore a matter of computing 
ai(X) and &2{X) with sufficiently enough bits. Given X = (Xi,...,X„) where Xi € {0,1} for 
i = l,...,n, the terms involve in ai{X) and dt2{X) are denoted by Cj = cos (uj — ^Xj^ and 
Sj = — sin (ujj — |Xj) and by Cj{k) and Sj{k) their k^^ their approximant obtained by truncating cj 
and Sj respectively to k bits. Since —1 < Cj < 1, —1 < ,Sj < 1, if we use 9 for possibly different real 
numbers at most one in absolute value, then we have cj = Cj{k) + sj = Sj{k) + The following 
algorithm computes troncated products &i{X) and a,2{X) in parallel. For the sake of clarity, in the 
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Algorithm 3 Computing the first m bits of tiie products C1C2 . . .Cn = ai(-X^) and S1S2 • • • = a2{X) 
in parallel using a communication tree configuration 

1: Every leaf of the tree, or node of level [log2(ra)], contains cj and Sj for j = 1, . . . ,n. 

2: Z •<— m + [log2(n)] (m will be given later) 

3: i ^ 1 
4: repeat 

5: We are at level [log2(ra)] — j- The are nodes, each containing a partial troncated product. 

In parallel, for i = 1, . . . , , person (i — 1)2^ + 1, located at node i of the current level, 
gets 2(1 — j) bits from the partial products {{I — j) bits for each product) that its right child 
has and two more bits for the signs of these two troncated products. 

In parallel, for i = 1, . . . ,n/2^ , person {i — 1)2^ + 1, located at node i of the current level 
computes the troncated products ^ = C(j_i)2J+i(^ ~ j)'^i2j{^ ~ j) ^ — ■^(i-i)2J+i(^ ~ 
J>j2jO - j)- 

The computed products Jl» and 4 are accurate to I — j — 1 bits. 

6: J J + 1 (go one level up) 
7: until j = flog2(ra)l. 



following description of the algorithm, assume the n is a power of two and the last level of the tree 
is full. 

By requiring person 1 to be located at the leftmost node of the tree, person 1 ends up being 
the chief at the root. At the end of the algorithm, person 1 knows exactly the first m bits of the 
products ai{X) and a.2{X). In fact, the original Cj's and Sj's being bounded by at most 1 in absolute 
value, a product of two paired troncated products is also bounded by 1 in absolute value; therefore 
1 bit is lost per computed product on a given level. The chief (root) knows (I — [log2(n)]) = m bits 
exactly from ai(X) and a2(X). The expected time of the algorithm in the bit complexity model is 
given by 

riog2(")i 

2(m + z + 1) = 2(m + l)riog2(n)l + [logs (n)l( [logs (n)l + 1). 

i=l 

Notice that the same algorithm can be used to simulate a Bernoulli random variable with 
parameter cos^(^^^) as given in the definition of the GHZ distribution. The n leaves are just 
angles {^pj + 4>j) + TjTr/2 for j' = 1, . . . , n, and we compute partial sums instead of partial products. 
Two bits arc expected from the parameters for generating a random Bernoulli variable so with the 
expected time using the previous algorithm is 0(log|(n)). 

The protocol 

Using the work of the previous section, we give an informal protocol for simulating the GHZ dis- 
tribution in the random bit model when all the parameters are spread out through the n persons, 
[protocol next page] 
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Algorithm 4 Protocol for sampling the GHZ distribution 
1: repeat 

2: The chief, person 1, generates a random bit B according to a Bernoulh random distribution 

with parameter cos^(^^^) using algorithm 3 suitably modified. 
3: if B = 1 then 

4: The chief adds vr to its w-parameter i.e. ui i— uji + ir (sampling P2). 
5: else 

6: The chief does not change its w-parameter (sampling pi). 
7: end if 

8: The chief generates a fair random bit S and broadcasts it to the other persons. 

9: Locally, everyone generates a random bit Xj according to a Bernoulli random distribution 

with parameter cos"^ (ujj). 
10: if 5" = 1 then 
11: Everyone flips their Xj. 
12: else 

13: They don't flip their Xj. 

14: end if 

15: Set j ^ 1. 

16: repeat forever 

17: Set k <(- j{n + 4) and m A; + 2. 

18: The chief computes the m-bits approximants of ai(X) and Si2{X) with algorithm 3 from 
the previous section and use them for computing the A;-bits approximants and Rk of 
U{qi{X) + q2{X)) and p{X) respectively. 

19: if Lfc -i?fc < then 

20: Set y ^ 1 and break. 

21: else a Lk — Rk > ^ then 

22: Set y ^ and break. 

23: else 

24: Set j ■(^ j + 1 and continue. 

25: end if 

26: end repeat 

27: until Y = 1 (accepting) 

28: The chief requires the other persons to broadcast their Xj and return X = {Xi, . . . 



The part contains within the "repeat forever" is essentially algorithm 2 from the previous section, 
and, as it was shown, it always stops. If J denotes the random number of iterations to produce 
y = or y = 1 and K denotes the number of bits needed to decide between acceptation (y = 1) and 
rejection (Y = 0), then E{J} = 1 since K = J(n + 4) and from the previous section EjX} = n + 4. 
The protocol uses algorithm 3 with an initial number of (n + 4) bits of precision for U{qi(X)+q2{X)) 
and p{X) that is the expected number of bits required to take a decision. 

Random bits are B, S and Xi,. . . ,X„ and they take an expected number of 2 + 1 + 2n bits. 
All other bits in the protocol are parameter-related bits. For producing B, an expected 0(log^(n) 
bits are communicated as explained from the previous section. The "repeat forever" requires an 
expected 0(nlog2(n)) bits of communication. The overall expected number of communicated bits 
is 0(nlog2(n)). 



9 



Discussion 



Knuth and Yao (1976) were the first to study the bit complexity for generating random integers 
with given distribution p{x) assuming a perfect BernoulU source of bits. They showed that any 
algorithm must use an expected number of bits at least equal to the entropy YlxPi-'^) lo§2(l/p(^))) 
and that the best algorithm does not use more than three extra bits above the lower bound. 

The GHZ distribution has binary entropy between and n. This provides the theoretical 
backdrop. Our problem, by virtue of the set-up in terms of parameters, whose bits must also 
sparingly be revealed, is slightly different, but related. 

For further results on the bit model in random variate generation, see chapter XIV of Devroye 
(1986) and the recent thesis of Lumbroso (2012). 
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